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The  electromagnetic  field  equations  and  Dirac  equations  for  oppositely  charged  wave 
functions  are  numerically  time-integrated  using  a  spatial  Fourier  method.  The  numerical 
approach  used,  a  spectral  transform  technique,  is  based  on  a  continuum  representation  of 
physical  space.  The  coupled  classical  field  equations  contain  a  dimensionless  parameter 
which  sets  the  strength  of  the  nonlinear  interaction  (as  the  parameter  increases,  interaction 
volume  decreases).  For  a  parameter  value  of  unity,  highly  nonlinear  behavior  in  the  time- 
evolution  of  an  individual  wave  function,  analogous  to  ideal  fluid  turbulence,  is  observed.  In 
the  truncated  Fourier  representation  which  is  numerically  implemented  here,  the  quantum 
turbulence  is  homogeneous  but  anisotropic  and  manifests  itself  in  the  nonlinear  evolution  of 
equilibrium  modal  spatial  spectra  for  the  probability  density  of  each  particle  and  also  for  the 
electromagnetic  energy  density.  The  results  show  that  nonlinearly  interacting  fermionic  wave 
functions  quickly  approach  a  multi-mode,  dynamic  equilibrium  state,  and  that  this  state  can  be 
determined  by  numerical  means. 
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1.  Introduction 

A  direct  numerical  simulation  of  the  self-consistent  electromagnetic  interaction 
between  two  oppositely  charged  and  densely  packed  spin  1/2  Dirac  particles  is  presented 
here.  While  the  approach  taken  here  is  nonperturbative,  it  is  not  based  on  lattice  gauge 
theory.  Instead,  it  is  a  solution  of  the  basic  set  of  coupled,  nonlinear  partial  differential 
equations  which  describe  a  fundamental  quantum  mechanical  system  in  terms  of  classical 
field  theory.  These  equations  are  integrated  forward  in  time  to  simulate  the  nonlinear 
evolution  of  two  Dirac  wave  functions  and  the  electromagnetic  field  which  couples  them. 

The  means  of  solution  is  a  Fourier  spectral  method,  in  which  a  three-dimensional 
momentum  space  (k-space)  is  restricted  to  contain  only  a  finite  number  of  discrete  modes 
(z.e.,  Ikl<k„„).  Although  k-space  is  discretized,  position  space  (x-space)  is  not:  The 
underlying  physical  space  is  a  continuum  and  not  a  discrete  set  of  points.  The  Fourier 
spectral  method  has  been  used  to  great  advantage  in  pioneering  work  in  turbulent  flow 
simulation  [1]  and  continues  to  be  used  in  the  study  of  turbulence  and  other  nonlinear 
dynamic  phenomena. 

Since  the  interaction  of  strongly  coupled  fields  is  essentially  nonlinear,  it  is 
generally  not  possible  to  assign  a  specific  frequency  to  each  spatial  mode.  The  actual  time 
dependence  of  the  Fourier  modes  is  found  by  integrating  the  equations  of  motion,  after 
which  a  time  sequence  for  each  mode  may  be  analyzed  to  determine  its  frequency  content 
during  the  sampling  time.  However,  it  must  be  remembered  that  in  the  evoultion  of  a 
nonlinear  dynamic  system  the  frequency  content  of  a  mode  is  constantly  changing,  as  the 
various  modes  are  nonlinearly  interacting  with  one  another.  To  robustly  determine  the 
frequency  spectrum  for  any  mode  requires  that  a  simulation  be  run  considerably  past  the 
time  at  which  the  initial  conditions  are  'forgotten'  by  the  nonlinear  system.  That  is 
computationally  expensive  (for  a  large  number  of  grid  points)  and  will  not  be  done  here, 
although  the  simulation  will  be  run  long  enough  to  see  the  establishment  of  spatial 
equilibria  for  the  interacting  wave  functions. 
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Strongly  interacting  quantum  mechanical  wave  fields  can  exhibit  the  same 
interesting  nonlinear  dynamic  behavior  seen  in  fluids  and  plasmas,  i.e.,  chaos  and 
turbulence.  (Chaos  and  turbulence  are  related  in  that  turbulence  may  be  thought  of  as  many 
degree-of-freedom  chaotic  motion,  while  "classical"  chaos  appears,  for  example,  when  the 
mathematical  model  of  a  turbulent  hydrodynamic  system  is  reduced  to  a  minimum  number 
of  degrees-of-freedom  [2].)  Here,  the  continuous  wave  functions  and  electromagnetic  field 
play  a  role  similar  to  that  of  conserved  components  in  a  mixture  of  classical  fluids;  for 
example,  probability  densities  are  analogous  to  component  mass  densities  as  both  satisfy 
identical  continuity  equations.  Turbulent  behavior  will  be  seen  in  the  dynamic  transfer  of 
energy  and  probability  between  different  spatial  modes  and  in  the  establishment  of  apparent 
equilibrium  modal  spatial  spectra. 

Again,  it  is  the  electromagnetic  interaction  of  oppositely  charged  spin- 1/2  particles 
(an  "electron"  and  a  "positron")  which  we  examine  here.  This  classical  "lepton-photon" 
system  is  described  by  Dirac  equations  and  the  electromagnetic  field  equations.  Although 
replacing  the  Dirac  equations  with  Schrodinger  equations  also  produces  a  system  which 
contains  the  nonlinear  electromagnetic  interaction,  the  coupling  parameter  is  small  and  the 
nonlinear  interaction  is  weak.  The  coupling  parameter  will  be  seen  to  increase  as  the 
density  of  particles  increases;  concomitantly,  the  mean  particle  velocity  will  become  more 
and  more  relativistic.  Thus,  the  description  of  a  lepton-photon  (or  any  fermion-gauge  field) 
system  with  a  strong  nonlinear  interaction  requires  the  use  of  the  Dirac  equation,  rather  than 
the  Schrodinger  equation. 

In  this  paper  a  first-quantized  or  "classical"  field  description  will  be  utilized,  which 
will  allow  us  to  follow  the  self-consistent  evolution  of  the  oppositely  charged,  two  particle 
quantum  mechanical  system.  The  basic  equations  will  be  given  in  nondimensional  form, 
followed  by  the  classical  Noether  invariants  of  the  system.  Then  the  numerical  method  will 
be  described,  and  numerical  results  will  be  presented. 
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In  addition  to  observing  turbulence  in  a  quantum  mechanical  system,  the  novelty  of 
the  work  presented  here  is  that  it  introduces  a  new,  nonperturpative  and  direct  approach  to 
studying  the  gauge  field  interaction  of  closely  packed  particles,  such  as  those  in  extremely 
dense  matter.  The  current  historical  context  is  similar  to  that  encountered  when  time- 
integration  methods  on  a  spatial  grid  were  introduced  into  the  study  of  general  relativistic 
flow  problems  [3]  and  into  nonrelativistic  quantum  processes  [4],  i.e.,  although  previous 
analytical  and  numerical  techniques  have  produced  and  continue  to  produce  many  valuable 
results,  time  integration  methods  allow  the  problem  at  hand  to  be  solved  (and  visualized) 
directly. 

A  study  of  coupled  nonlinear  Dirac  equations  in  four  dimensions  has  appeared 
before,  in  the  work  of  Alvarez  [5],  where  soliton-like  behavior  was  examined.  In  that 
work,  however,  the  mediating  gauge  field  was  eliminated  by  introducing  ad  hoc  terms  into 
two  separate  free-particle  Dirac  equations  so  as  to  produce  a  direct  nonlinear  coupling. 
Here,  we  study  two  classical  Dirac  fields  realistically  coupled  by  an  electromagnetic  field 
and,  in  this  case,  do  not  find  the  'blow-up'  problem  which  appears  in  the  direct  nonlinear 
coupling  model  [6].  The  approach  taken  here  is  also  generically  similary  to  that  of 
Bialynicki-Birula,  et  al.,  who  have  recently  examined  the  self-consistent  time  evolution  of 
quantum  fields  in  terms  of  the  Wigner  distribution  function  [7]. 
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2.  Basic  Equations 


The  Dirac  equation  is  (here,  standard  notation  [8]  is  used) 

JT  =  mcS-,  p,  ^ 


(1) 


For  electrons,  e  is  the  negative  electronic  charge  and  for  positrons,  e  is  the  positive 
electronic  charge;  m  is  the  electron  or  positron  mass,  c  is  the  speed  of  light,  and  h  is 
Planck’s  constant.  Explicitly,  the  (4x4)  Dirac  matrices  are 


yO  =  P,  Y=pa, 


(2) 


(Greek  indices  range  from  0  to  3  with  a  metric  signature  of  +— ,  while  Latin  indices  range 
over  1,  2,  3  (/.  e.,  x,  y,  z)  with  a  metric  signature  of  +++;  repeated  indices  imply 
summation.  Also,  boldface  denotes  a  3-vector.) 

The  electron  and  positron  wave  functions  are  complex  entities  and  will  be  expressed 

here  as 
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where  the  /?/,  S/,  Pi,  and  Qi  (/=!, 2,3,4)  are  real  functions  of  time  and  position.  Coupled 
with  the  Dirac  wave  equations  are  the  electromagnetic  field  equations  (using  the  Lorentz 
gauge  condition): 
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= 0 


jS  =  nY'‘'5'p.  Te>p'V 

The  conservation  of  probability  is  guaranteed  by  the  continuity  equations  each  wave 
function  satisfies:  ^^j^‘=^lP+V•j=0,  where  =  (p  j)  corresponds  to  either  particle. 

At  this  point,  we  will  nondimensionalize  equations  (1)  and  (4).  Since  a  spatial 
Fourier  method  will  presently  be  used  for  numerical  simulation,  units  of  distance  will  be  in 
terms  of  Lo/27t,  where  is  the  side  length  of  the  periodic  physical  cube.  Using  the  operator 
equivalence  given  in  (1),  the  nonlinearly  coupled,  nondimensional  dynamic  equations  are 

•yt'lap+iAjT,  =  -i'P.,  f  (ap-iApj'Pp  =  -iTp 

(af-vV^-KCiS-j?).  apAi*=o  ,5) 


These  equations  contain  only  one  parameter  which  determines  the  nature  of  the  interaction: 


where  the  Compton  wavelength  of  the  electron  is  Xc  =  h/mc  =  2.43  pm  and  the  fine-structure 
constant  is  a  =  27te^/hc  =  1/137.  Here,  'F*  (s=e,p)  is  normalized  so  that  the  integral  of  js 
over  the  characteristic  volume  L<,3  is  equal  to  unity  Note  that  for  the  interaction  parameter  to 
have  a  value  of  unity  (k  =  1),  then  Lo  =  =  0.47  pm  and  the  density  of  particles  must  be 

around  Lo'^  =  10^’  cm-^;  electrons  at  densities  up  to  10^’  cm  ^  are  believed  to  exist  in  the 
outer  layers  of  neutron  stars  [9].  This  density  is  also  achieved  by  scattering  particles  whose 
'interaction  time'  is  at  least  10  ^'  seconds,  i.e.,  'resonant'  particles. 
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3.  Noether  Invariants 

The  classical  invariants  [10]  of  the  electromagnetically  interacting  electron-positron 
system  can  be  derived  from  the  Lagrangian  density 


A  =  +  TpY'‘{l^'J'p)  -  (OpYplyf  >Pp] 

-  T.y.  -  Tp-Fp  -  ^FpvFK''  (7, 

where  the  "covariant"  derivative  Du  and  the  electromagnetic  field  tensor  are, 
respectively, 

=  3|^+iA^,  Fjiv  ^  3jxAy-3v  Ajj^. 

The  nondimensional  volume  of  the  3-space  cube  is  (27t)^;  an  integral  over  this  volume  is 

(Q)  =  l^Q(t,r)dr 

a  definition  which  allows  for  notational  conciseness.  Using  Noether's  theory  [10]  (along 
with  K  =  1),  the  important  classical  invariants  are  found  to  be 

Normalization  (total  probability): 
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Energy: 


E  =  (-i'FeY-  V  -i^pY-  V  'Fp  +^p^p  ‘A  -Oe-jp)  +  ^(l  E|2+|  H \^)) 

(E^-atA-VA°,  H  =  VxA,  je^^eY%,  jp^^^pY^p) 


(9) 


Momentum: 


P  =  (p)  =  (-i^eY°  V  ‘Fe  -i^pf  V  'Fp  +A{  j°  -  ji)  +  Exh) 


(10) 


Angular  Momentum: 


=(' 


rxp  +^4^6  Y°  2:  'Fe  +^'FpY°2:H^p  +  Ax 


(11) 


The  invariants  (8)  through  (1 1)  are  classical  and  the  fields  contained  in  them  are  not 
considered  to  consist  of  explicit  creation  and  annihilation  operators  [10];  according  to  the 
tenets  of  Lagrangian  field  theory,  these  invariants  should  be  preserved  during  the  time- 
evolution  of  a  closed  system.  The  invariance  of  a  numerical  model  based  on  equations  (5) 
will  be  examined  in  the  next  section.  (The  specific  parts  of  these  invariants  which  are 
associated  with  either  the  electron,  positron,  or  photon  fields,  or  with  their  interaction,  can 
be  easily  separated  out  of  the  total  expression  and  examined  individually,  as  required.) 
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4.  Numerical  Results 

Using  the  equations  given  in  (5),  the  time-evolution  of  the  lepton-photon  system 
was  simulated  on  a  32x32x32  k-space  grid.  Simulations  were  performed  using  a  spatial 
Fourier  transform  method  [11]  with  a  de-aliased  [12],  third-order  time-integration  scheme 
[13]  (this  approach  is  somewhat  similar  to  that  used  in  nonrelativistic  quantum  mechanics 
simulations  [14]).  For  numerical  simulation,  the  equations  (5)  are  expressed  as  follows: 

=  -[a-V  i(p  (p  -  a-A)]'Fe,  at'I'p  =  -[ol-V  -h  i((3  -  tp  -h  a  AjJ'Fp 

a^c  =  v'a K(  je  -  jp),  atA  =  c,  at(p  =  -VA 

(P  =  A°,  jp  ^ 'I^pY'I^p  (12) 

If  we  set  K=0  and  assume  that  9  and  A  are  initially  zero,  then  the  equations  for  'Fe  and 
are  linear.  In  this  case,  both  'F*  and  %  have  linear  solutions  ('F  being  either  one): 

'F(x,t)  =  ^  [cos(cOkt)  -  i(ojj^(p  +a  k)sin(tOkt)]o(k)  e'*^  * 

k  (13) 

where  0(k)  is  a  time-independent,  complex,  four  component  column  vector  and  to^  = 
(k^+l)'^.  The  lowest  frequency  is  obviously  cOo  =  1,  with  a  corresponding  period  of  = 
27t.  Even  though  we  will  examine  non-linear  behavior  (k>0)  and  will  not  utilize  (13) 
further,  (13)  indicates  that  a  simulation  needs  to  be  run  from  t=0  to  at  least  t=To. 

The  fields  which  comprise  the  system  (12)  are  seen  to  be,  using  (3), 

Ti,  R2,  R3,  R4,  Si,  S2,  S3,  S4 

Ax,  Ay,  Aj,  Cx>  Cy,  Cz 

P].  Pi,  P3,  P4,  Qi,  Qi,  Q3.  Q4  (14) 
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In  the  numerical  method,  these  are  expanded  in  terms  of  spatial  Fourier  series,  for  example. 


Ri(x,t)  =  XRi(k,t)  eik  * 

k  (15) 

Thus,  the  few  non-linear  partial  differential  equations  in  (12)  are  transformed  into  many 
non-linear  ordinary  (in  time)  differential  equations. 

In  addition  to  the  equations  in  (12),  there  is  also  an  auxiliary  condition  which  must 
be  satisfied: 


-V^(p  =  K(pe  -  pp)-t-V  C 

where 

=  Pp  =  = 'FpP'Fp 

Equation  (16)  arises  when  the  Lorentz  condition  and  the  wave  equation  for  the  electric 
potential  <p  are  combined.  Thus,  in  (14)  (p  is  not  an  independent  dynamic  variable; 
however,  either  (16)  or  the  Lorentz  gauge  dt(p  =  -V  -  A  can  be  used  to  determine  (p  during 
the  dynamic  evolution  of  the  system  (whichever  is  more  computationally  efficient  -  here  the 
Lorentz  condition  is  used).  Initially,  however,  (16)  is  always  needed  to  determine  (p. 

In  a  spatial  Fourier  method,  the  Lorentz  condition  is  d(p(k)/dt  +ik  - A  (k)  =  0,  and 
a  gauge  transformation  of  the  electromagnetic  fields  has  the  modal  form  {A(k),  C(k)}— > 
{ A(k)  -  ik0(k),  C(k)  -  ikd9(k)/dt}.  Here  0(k)  satisfies  the  modal  wave  equation 
d20(k)/dt2  +  k20(k)  =  0.  Under  a  gauge  transformation,  the  modal  form  of  the  change  of 
the  quantum  mecahnical  wave  function  is  'F(k)->  {exp(i0)4^}(k);  i.e.,  the  modal  gauge 
transformation  is  just  the  spatial  Fourier  transform  of  the  physical  space  gauge 
transformation,  as  long  as  all  possible  modes  k  are  retained.  This  last  stipulation  results 
from  the  observation  that  the  spatial  Fourier  transform  of  exp(i0)  must  have  an  infinite 
number  of  modes,  and  the  numerical  method  cannot  de-alias  a  quadratic  product  where  one 


(16) 

(17) 
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of  the  cofactors  is  known  to  contain  more  than  the  truncated  set  of  modes.  However,  if 
invariance  under  only  an  infinitesimal  gauge  transformation  'F(k)^  {(1+  i0)'P}(k)  is 
required  (at  each  numerical  time-step),  and  0  is  restricted  to  contain  no  modes  outside  the 
truncated  set,  then  the  numerical  method  is  gauge  invariant. 

In  the  present  numerical  method,  each  field  (for  an  NxNxN  spatial  grid)  has 
approximately  0.4388  degrees-of-freedom  (real  and  imaginary  parts  of  independent 
Fourier  modes).  Thus  each  of  the  fields  in  (14)  for  a  32^  grid  has  about  14400  degrees-of- 
freedom,  and  since  there  are  22  independent  fields,  the  model  system  has  a  total  of  about 
320,000  degrees-of-freedom.  Computationally,  this  means  that  we  have  this  many 
coupled,  nonlinear,  ordinary  differential  equations  to  solve.  (Although  it  will  not  be  done 
here,  an  an  N=64  simulation  is  also  be  possible.  This,  of  course,  requires  a 
correspondingly  greater  investment  of  computer  resources.) 

One  long  simulation  will  be  presented  in  detail  here.  This  was  run  on  a  Cray  YMP, 
with  a  cpu  time  per  simulation  time  step  of  about  14  seconds  for  N=32.  The  coupling 
constant  was  k  =  1  and  the  initial  conditions  were  such  that  only  R]  S2,  P2,  and  Q,  were 
nonzero  with  <Ri2>=<S2^>=<  P22>=<Qi2>  {i.e.,  neither  the  electron  nor  the  positron  had 
any  initial  kinetic  energy,  linear  momentum,  or  angular  momentum).  Initially,  Rj  S2,  P2, 
and  Qi  were  described  by  spatial  three-dimensional  gaussian  density  distributions  centered 
on  the  grid  points  (8,16,16),  (16,8,8),  (24,16,16),  and  (16,24,16),  respectively;  all  had 
standard  deviations  of  8  grid  spacings  and  were  set  to  zero  beyond  one  standard  deviation 
from  their  respective  centers.  Also,  at  t  =  0,  A  =  C  =  0,  and  (|>  was  determined  by  (16). 
Each  computational  time  step  advanced  the  system  At=0.{)()0125  simulation  time  units. 

During  the  simulation,  which  ran  from  t=  0  to  6.3  {i.e.,  2n),  the  normalization  of 
the  electron  and  positron  wave  functions  was  conserved  to  1  part  in  10*  (thus,  total  charge 
was  conserved  to  this  accuracy,  and  there  was  no  'blow-up'  [6]).  The  total  energy  given  in 
(9)  was  also  conserved  extremely  well,  fluctuating  no  more  than  0.04  %  during  the  run. 


Thus,  the  Noether  invariants  of  nomalization  (/.  e.,  total  charge,  probability,  or  particle 
number)  and  energy  were  essentially  conserved  during  the  run. 

Another  measure  of  numerical  efficacy  lies  in  behavior  of  the  "center  of  inertia"  R 
of  the  system,  which  should  remain  fixed  (since  the  P=0  at  t=0)  for  both  runs.  Here,  v/e 
define  R  as 


where  P  is  the  total  momentum,  as  given  by  (10),  and  E  is  the  total  energy,  as  given  by 
(9).  Since  the  edge  length  of  the  computational  3-D  volume  is  2k,  the  percent  variation  is 
defined  as  100%xlR(t)l/27c.  The  fluctuation  in  the  center  of  inertia  was  less  than  0.4  %, 
commensurate  with  the  numerical  variation  in  energy. 

To  get  an  appreciation  of  the  difference  between  linear  and  non-linear  evolution, 
consider  Fiqure  1,  where  the  linear  and  non-linear  time  dependence  of  the  Fourier 
coefficients  Ri(k,t)  and  Si(k,t)  for  k=0  is  compared  (for  k=0,  all  coefficients  are  real). 
According  to  (13)  the  linear  behavior  of  the  pair  should  be  Ri(0,t)=cos(t)  and  Si(0,t)=- 
sin(t)  (for  this  figure,  the  amplitude  has  been  normalized  to  unity).  The  actual  trajectory  of 
the  pair  is  obviously  different  from  the  linear  prediction;  there  are  clearly  many  more 
frequencies  present  than  just  the  single  one  corresponding  to  the  linear  mode.  In  fact,  if  the 
behavior  of  any  coefficient  is  examined  in  a  similar  manner,  the  same  behavior  will  be  seen: 
a  'random  walk'  around  the  origin. 

To  get  another  view  on  the  dynamic  evolution  of  the  model  system,  let  us  break  up 
the  total  energy  (9)  into  its  constituent  parts: 

Ee  =(-i'PeY- V  'Pe  Ep  =  (- i^ p Y’ V  'Fp  -»-'Pp'Pp) 

£l=(-A-(jc-jp)),  EEM  =  (^(|E|2-h|H|2)) 

1  1 


(19) 


Here  we  have  defined  the  "electron's  energy"  Eq,  "positron's  energy"  E^,  "interaction 
energy"  Ei,  and  electromagnetic  energy  E^^.  The  evolution  of  these  energies  is  shown  in 
Figure  2. 


Next,  consider  the  quantities 


(|'Fe.iP)  =  (Ri  +S?),  (|'Fp.ip)  =  (p?  +Qf),  i  =  1,2, 3,4 


(20) 


These  are  just  the  contributions  each  component  makes  to  their  respective  normalization 
integral  (8).  Their  time  evolution  is  given  in  Figures  3  and  4  for  the  electron  and  positron 
fields,  respectively. 


5.  Quantum  Mechanical  Turbulence 

Let  us  now  take  up  the  matter  of  nonlinear  dynamics  and  turbulence  in  this  multi- 
mode  quantum  mechanical  system.  The  parameter  k  plays  a  role  in  (12)  analogous  to  that 
played  by  the  Reynolds  number  in  fluid  turbulence,  i.  e.,  as  these  numbers  tend  to  zero  the 
nonlinear  effects  in  the  respective  systems  disappear.  A  characteristic  of  turbulent  behavior 
is  the  manner  in  which  energy  (or  probability)  is  shared  nonlinearly  between  different 
modes.  This  is  illustrated  for  the  present  simulation  in  Figures  5  and  6,  where  the  wave 
number  spectra  of  the  electron  probability  density  and  electromagnetic  energy  density, 
respectively,  at  several  different  times  are  shown.  (The  positron  spectra  are  very  similar  to 
the  electron  spectra.)  (In  a  32^  run,  the  maximum  wave  number  is  15.07,  and  thus 
log(k„,„)=1.17).  The  spectra  shown  in  Figures  5  and  6  are  derived  from  the  modal 
densities  by  finding  the  average  over  all  k  with  the  same  magnitude  lkl=k,  and  multiplying 
this  average  by  k^.  These  "omnidirectional"  spectra  are,  explicitly: 

ji[R?(k)+S?(k)]| 

^2(2^  k^Af(k)  -KC?(k)]| 

P(k)  =  i;|^S  (t  [P?(k)+Q?(k)]| 

I  (21) 

Here,  N(k)  is  the  number  of  terms  in  the  summation  over  k  such  that  lkl=k.  (The  values  of 
Pe(k)  and  Pem(k)  shown  in  Figures  5  and  6,  respectively,  are  smoothed  by  averaging  over 
nearest  neighbors.) 

The  shape  of  the  spectra  at  t  =  0  is  the  initial  spectra  (in  a  linear  run,  where  k=0, 
these  spectra  do  not  change  shape  at  all  with  time).  As  is  seen  in  Figures  5  and  6,  there 
was  a  considerable  amount  of  energy  and  probability  transfer  between  the  different  modes; 
in  fact,  all  the  spectra  appear  to  be  converging  to  equilibria. 
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It  should  be  possible  to  predict  these  spectra  a  priori  as  is  done,  for  example,  in 
ideal  three-dimensional  magneto- fluid  turbulence  [15],  since  the  system  of  equations  (12) 
satisfy  all  the  criteria  necessary  for  'absoulte  equibrium  ensemble'  theory  to  apply  [16].  In 
particular,  a  partition  function  involving  the  numerical  invariants  of  the  model  system  is 
determined  and  used  to  construct  canonical  ensemble  predictions,  for  example,  of  turbulent 
energy  spectra  [17].  (This  procedure  has  a  close  analogy  to  work  in  lattice  field  theory, 
where  a  partition  function  involving  a  Eucidean  action  is  sought  [18];  remember  though, 
that  the  model  underlying  the  simulation  here  is  a  continuum  model,  while  that  of  lattice 
gauge  theory  is  not.)  However,  in  non-dissipative  fluid  turbulence,  the  invariants  are 
quadratic  sums,  while  the  situation  here  is  more  complicated  since,  for  example,  the 
interaction  energy  Ei  is  cubic  in  nature,  and  the  relation  (16)  between  the  potential  and  the 
dynamical  fields  introduces  a  term  quartic  in  the  wave  functions  into  the  energy  expression. 
Developing  this  possibility  will  be  deferred. 

In  order  to  actually  "see"  the  interaction,  consider  Figure  7.  This  figure  indicates 
the  relative  values  of  the  electron's  3-D  probability  density  the  electromagnetic  3- 

D  energy  density  femCx).  and  the  positron's  3-D  probability  density  l'Pp(x)P,  summed  in 
the  z-direction  and  projected  onto  the  x-y  plane  for  equally  spaced  times  during  the  run. 

Figures  5,  6,  and  7  indicate  a  relative  change  in  size  of  the  various  physical  fields 
with  time,  a  change  which  can  be  quantified  by  defining  wave  numbers  Kep  and  Kem  : 


X  k2[/'e(k)-^Pp(k)]  X  k^/^emCk) 

K2  _  k  vZ  _  k 

X  [/’e(k)+Fp(k)]  X  ^em(k) 

k  k  (22) 

where  the  (i=e,p,em)  are  given  in  (21).  The  time  evolution  of  these  root-mean-square 
(rms)  wave  numbers  are  shown  in  Figure  8. 


Although  the  spectra  shown  in  Figures  5  and  6  and  defined  by  (21),  and  the  rms 
wave  numbers  shown  in  Figure  8  and  defined  by  (22),  are  determined  by  averaging  over 


all  directions,  the  turbulence  which  is  simulated  is  not  in  fact  isotropic.  We  can  define  a 
measure  of  anisotropy  in  the  x-direction  as  follows: 

Then  measures  of  anisotropy  in  the  y-  and  z-directions  are  My=Nyzx  and  Mz=Nzxy, 
respectively,  and  satisfy  Mx+My+Mz=0  (these  measures  are  similar  to  those  used  in  fluid 
turbulence  work  [19]).  The  quantity  'F  can  be  either  the  electron  or  positron  wave  function 
or  the  complex  electromagnetic  vector  A+iC.  The  quantities  Mx,  My,  and  Mz  change  with 
time;  in  Figure  9,  the  evolution  of  these  quantities  for  the  electron  (4'='Fe)  is  shown. 
Measures  of  the  positron  and  electromagnetic  anisotropy  behaved  very  similarly. 

The  observed  anisotropy  occurs  because  we  have  a  mixture  of  "charged  fluids".  At 
t=0  the  electron  and  positron  densities  are  separated,  more  or  less,  along  the  x-axis  and  are 
initially  motionless.  Since  the  initial  densities  are  composed  of  spherical  distributions,  and 
have  no  motion,  we  have  Mx=My=Mz=0  at  t=0,  according  to  (23),  for  the  elctron  and 
positron  (the  electromagnetic  field  has  a  slight  initial  anisotropy).  However,  as  the  particles 
are  electrically  attracted,  they  begin  to  "move"  in  response  to  one  another,  and  this  is 
reflected  by  gradients  in  the  x-direction  increasing  more  quickly  than  gradients  in  the  other 
two  directions.  Hence,  the  anisotropy  is  contained  in  the  initial  conditions  and  manifests 
itself  in  the  direction  of  "plasma  oscillations".  Thus  we  have  homogeneous  (because  the 
nonlinear  dynamics  occurs  in  a  periodic  cell  in  space)  but  anisotropic  turbulent  phenomena. 


6.  Conclusion 

In  this  article,  the  numerical  simulation  of  a  nonlinear  quantum  mechanical  lepton- 
photon  interaction  has  been  described.  This  simulation  was  performed  using  a  truncated 
spatial  Fourier  method  to  march  the  classical  system  equations  forward  in  time.  It  was  seen 
that  the  interaction  was  highly  nonlinear,  and  that  the  time-evolution  of  the  wave  functions 
could  be  described  as  'turbulent',  when  the  interaction  was  nominally  "strong"  (t.e.,  k=1). 
Conversely,  if  setting  k=1  results  in  the  rapid  transfer  of  probability  and  energy  between 
spatial  modes,  then  this  sets  the  natural  'equilibrium'  interaction  scale  length  as 
(=  0.47  pm  for  electrons).  That  these  closely  interacting  fermionic  wave  functions  quickly 
approached  a  multi-mode,  dynamic  equilibrium,  is  indicated  by  the  various  Figures. 

Natural  extensions  of  the  present  work  are  the  following.  First,  propagating 
electromagnetic  fields  can  be  introduced  into  the  initial  conditions  to  examine  their  effects 
on  the  behavior  of  the  system.  Second  higher  resolution  (e.  g.,  64^)  runs  can  be 
performed,  perhaps  on  a  massively  parallel  processing  system.  Third,  a  statistical  theory 
based  on  a  classical  partition  function  involving  system  invariants  can  be  developed. 
Fourth,  the  work  can  be  extended  to  encompass  nonabelian  gauge  fields.  This  last 
possibility  is  an  intriguing  one  as  it  could  provide  a  non-perturbative  method  for  studying 
few-body  interactions  in  such  quantum  systems  as  the  quark-gluon  plasma. 

The  classical  results  described  here  pertain  only  to  a  quantum  mechanical  system  of 
interacting  single  particles.  A  multiparticle  treatment  must,  of  course,  be  based  on  quantum 
field  theory,  as  has  been  done,  for  example,  by  Bialynicki-Birula,  et  al.  [7].  The  main 
intent  here,  however,  was  to  demonstrate  highly  nonlinear  behavior  in  a  quantum 
mechanical  system  and  to  present  a  numerical  method  for  observing  that  behavior.  In  so 
doing,  we  have  shown  that  a  microscopic  quantum  mechanical  system  and  macroscopic 
classical  system  (such  as  a  fluid)  have  a  common  mechanism,  i.e.,  a  parameteric  non-linear 
coupling,  which  can  induce  a  host  of  interesting  phenomena  (such  as  turbulence)  into  the 
dynamical  behavior  of  either  system. 
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Figure  2.  Time  variation  of  electron,  positron,  electromagnetic,  and  interaction  energies. 
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Figure  4.  Positron  components  i=l,2,3,4)  versus  time. 
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Figure  7.  Spatial  densities  of  electron  and  positron  probability,  and  electromagnetic  energy, 
averaged  over  z  and  projected  onto  x-y  plane,  at  various  times  (black-lowest  to  white-highest). 
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Figure  8.  RMS  wave  numbers  for  the  combined  electron-positron  spectra  (Kcp)  and 
for  the  electromagnetic  spectra  (Kem). 
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Figure  9.  Measures  of  anisotropy  for  tiie  turbulent  evolution  of  the  electron  wave  function.. 
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